//- update saturation 
theta = pcModel->correctAndSb(h);

//- relative permeability computation
krModel->correct();
krthetaf = fvc::interpolate(krtheta,"krtheta");

//- mobility computation 
Mf = rhotheta*mag(g)*Kf*krthetaf/mutheta;
Lf = rhotheta*Kf*krthetaf/mutheta;

//- compute fluxes
phiG = (Lf * g) & mesh.Sf();
phi = phiG-(Mf*fvc::snGrad(h))*mesh.magSf();
Utheta = fvc::reconstruct(phi);
